Read in real temp data and test function

Read in other temperatures and compare

#setwd("~/VegParamsSA/R/")
low_tmin = read.csv("../clim/temp_low.tmin")
low_tmax = read.csv("../clim/temp_low.tmax")
low = cbind(low_tmin, low_tmax)
colnames(low) = c("tmin","tmax")
low$date = seq.Date(from=as.Date("2009/10/1"), length.out=nrow(low), by="day")
low$tavg = (low$tmin+low$tmax)/2
low$m_tavg = leaf_cond_temp_curve(low$tavg)[[2]]

hot_tmin = read.csv("../clim/temp_extreme.tmin")
hot_tmax = read.csv("../clim/temp_extreme.tmax")
hot = cbind(hot_tmin, hot_tmax)
colnames(hot) = c("tmin","tmax")
hot$date = seq.Date(from=as.Date("2009/10/1"), length.out=nrow(hot), by="day")
hot$tavg = (hot$tmin+hot$tmax)/2
hot$m_tavg = leaf_cond_temp_curve(hot$tavg)[[2]]

ggplot(m_tavg_df)+geom_point(aes(x=tair, y=m_tavg2)) 

ggplot() + geom_point(data=low, aes(x=tavg, y=m_tavg, col='low')) + ggtitle("lower temps by 4.16 deg C")

ggplot() + geom_point(data=hot, aes(x=tavg, y=m_tavg, col='high')) + ggtitle("higher temps by 4.6 deg C")

ggplot(m_tavg_df)+
  geom_density(aes(x=tair, col='hist')) +
  geom_density(data=low, aes(x=tavg, col='low')) +
  geom_density(data=hot, aes(x=tavg, col='high')) +
  geom_vline(aes(xintercept=15, linetype='topt'))

m_tavg_df$hist = m_tavg_df$m_tavg2
m_tavg_df$low = low$m_tavg[low$date < as.Date("2018-12-25")]
m_tavg_df$hot = hot$m_tavg[hot$date < as.Date("2018-12-25")]

Look at scenario output of trans and psn

compare with the daytime temp leaf conductance multiplier

tmp = readout %>% dplyr::filter(code=="quag" & climmod !="hot") 
tmp$climmod[tmp$climmod=="hotter"] <- 'hot'
tmp$climmod[tmp$climmod=="lower"] <- 'low'

m_tavg_long <- pivot_longer(m_tavg_df, cols=c('hist','hot','low'), names_to="climmod", values_to="m_tavg") %>% 
  dplyr::select(date, climmod, m_tavg)

allout = inner_join(tmp, m_tavg_long, by=c("date","climmod"))

ggplot(allout) + geom_line(aes(x=date, y=m_tavg, col=climmod))

allout %>% group_by(yd, climmod) %>% summarize_at(vars(m_tavg), list(mean)) %>% ggplot() + geom_line(aes(x=yd, y=m_tavg, col=climmod)) + ggtitle("avg. daytime temp multiplier on leaf conductance")

ggplot(allout) + geom_point(aes(x=m_tavg, y=gpsn)) +
  geom_smooth(method='lm', aes(x=m_tavg, y=gpsn)) +
  facet_wrap('climmod')
## `geom_smooth()` using formula 'y ~ x'

ggplot(allout) + geom_point(aes(x=m_tavg, y=trans)) +
  geom_smooth(method='lm', aes(x=m_tavg, y=trans)) +
  facet_wrap('climmod')
## `geom_smooth()` using formula 'y ~ x'

ggplot(allout) + geom_density2d(aes(x=m_tavg, y=trans, col=climmod)) + facet_wrap('climmod')
## Warning: stat_contour(): Zero contours were generated
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf

ggplot(allout[allout$climmod=='hot',]) + geom_density2d(aes(x=m_tavg, y=trans)) + ggtitle("daily transpiration [mm] - hot")

ggplot(allout[allout$climmod=='hot',]) + geom_density2d(aes(x=m_tavg, y=gpsn)) + ggtitle("daily gpsn - hot")

topt_range = seq(from=12, to=22, by=1)
tcoeff_range = rnorm(100, mean=0.2, sd=(.2*.2))


m_tavg_topt = map_dfc(topt_range, ~leaf_cond_temp_curve(tair_avg$tair, topt=.x)[[2]])
colnames(m_tavg_topt) = topt_range
m_tavg_topt$tair = tair_avg$tair
m_tavg_topt_df = pivot_longer(m_tavg_topt, cols = 1:10, names_to="topt", values_to="m_tavg")
ggplot(m_tavg_topt_df) + geom_line(aes(x=tair, y=m_tavg, col=topt))

ggplot(m_tavg_topt_df) + geom_boxplot(aes(x=topt, y=m_tavg, col=topt))

m_tavg_tcoeff = map(tcoeff_range, ~leaf_cond_temp_curve( tair=tair_avg$tair, tcoef = .x)[[2]])

tcoef_df = as.data.frame(sapply(m_tavg_tcoeff, c))
colnames(tcoef_df) <- tcoeff_range
tcoef_df$tair = tair_avg$tair
tcoef_df2 = pivot_longer(as.data.frame(tcoef_df), cols = 1:100, names_to="tcoef", values_to="m_tavg")
tcoef_df2$tcoef = as.double(tcoef_df2$tcoef)

ggplot(tcoef_df2) + geom_point(aes(x=tair, y=m_tavg, col=tcoef)) +
  geom_line(data=tcoef_df2[tcoef_df2$tcoef>=0.1999 & tcoef_df2$tcoef<=0.201,], 
            aes(x=tair, y=m_tavg, col='default value 0.2'), col='red')

ggplot(tcoef_df2) + geom_point(aes(x=tcoef, y=m_tavg))

# LHS copied from ESM232
factors = c("tcoef","topt")
# Decide How many parameter sets to run
nsets=100
# choose distributions for parameters - this would come from
# what you know about the likely range of variation
q = c("qnorm", "qunif")
q.arg = list(list(mean=0.2,sd=0.2), list(min=12, max=22))

# generate samples from LHS
sens_leaf = LHS(NULL,factors,nsets,q,q.arg)
sens_pars = get.data(sens_leaf)
res = sens_pars %>% pmap(leaf_cond_temp_curve, tair=tair_avg$tair)
resd = as.data.frame(sapply(res,c)[1,])
sens_leaf2 = pse::tell(sens_leaf, as.matrix(resd), res.names="m_tavg")
pse::plotscatter(sens_leaf2, col="blue", cex=5)

# prcc's automatically generated and easy to plot
pse::plotprcc(sens_leaf2)

sens_leaf2$prcc
## [[1]]
## 
## Call:
## pcc.default(X = L, y = r, rank = T, nboot = nboot)
## 
## Partial Rank Correlation Coefficients (PRCC):
##         original
## tcoef -0.9337254
## topt   0.2746330